// 版权所有2010 Go作者。保留所有权利。
// 此源代码的使用受BSD样式
// 许可证的约束，该许可证可以在许可证文件中找到。

package math

/*
	Bessel function of the first and second kinds of order zero.
*/

// 原始C代码和下面的长注释是
// 来自FreeBSD的/usr/src/lib/msun/src/e_j0。c和
// 随此通知一同发出。go代码是原始C.
// 
// ==============================================================================================================================
// Sun Microsystems，Inc.版权所有。
// 
// 由Sun Microsystems，Inc.的企业SunPro开发。
// 使用、复制、修改和分发本
// 软件的权限是免费授予的，前提是保留本通知
// 的权利。
// /================================================================================================================================================
// /
// /第一类和第二类零阶贝塞尔函数。
// Method--j0（x）：
// 1。对于微小的x，我们使用j0（x）=1-x**2/4+x**4/64-。。。
// 2。将x减少到| x |，因为j0（x）=j0（-x），并且对于x in（0,2）
// j0（x）=1-z/4+z**2*R0/S0，其中z=x*x；
// （精度：|j0-1+z/4-z**2R0/S0 |<2**-63.67）
// /对于x in（2，inf）
// j0（x）=sqrt（2/（pi*x））*（p0（x）*cos（x0）-q0（x）*sin。最好计算sin（x0），cos（x0）=cos（x）cos（pi/4）+sin（x）sin（pi/4）
// /=1/sqrt（2）*（cos（x）+sin（x））
// sin x0）=sin sin（x）cos（pi/4）-cos（x）sin（pi/4）
// sin（x）+-cos（x）=-cos（2x）/（sin（x）-+cos（x））
// 计算更差的一个。）
// 
// 3个特例
// j0（nan）=nan 
// j0（0）=1 
// j0（inf）=0 
// 
// 方法--y0（x）：
// 1。对于x<2。自
// 其中
// V（z）=1+v01*z+…+v04*z**4 
// 绝对近似误差范围为2**-72。注意：对于微小的x，U/V=u0和j0（x）~1，因此
// y0（微小）=u0+（2/pi）*ln（微小），（选择微小<2**-27）
// 2。对于x>=2。
// y0（x）=sqrt（2/（pi*x））*（p0（x）*cos（x0）+q0（x）*sin（x0））
// 其中x0=x-pi/4。最好用上述方法计算sin（x0），cos（x0）
// 。
// 3。特殊情况：y0（0）=-inf，y0（x<0）=NaN，y0（inf）=0。
// 

// J0返回第一类零阶贝塞尔函数。
// 
// 特殊情况为：
// J0（±Inf）=0 
// J0（0）=1 
// J0（NaN）=NaN 
func J0(x float64) float64 {
	const (
		Huge   = 1e300
		TwoM27 = 1.0 / (1 << 27) // 2**-270x3e40000000000000000000000
		TwoM13 = 1.0 / (1 << 13) // 2**-130x3f200000000000000000
		Two129 = 1 << 129        // /2**1290x480000000
		R04 = 1.82954049532700665670e-06  // 0x3EBEBB1D10C503919 
		R05 = -4.61832688532103189199e-09 // 0xBE33D5E773D63FCE 
		S01 = 1.56191029464890010492e-02  // 0x3F8FFCE882C8C2A4 
		S02 = 1.16926784663337450260e-04  // 0x3F1EA6D2D57DBF4 
		S03 = 5.13546550207318111446e-07  // 0x3EA13B54CE84D5A9 
		S04 = 1.16614003333790000205e-09  // 0x3E1408BCF4745D8F 
	)
	switch {
	case IsNaN(x):
		return x
	case IsInf(x, 0):
		return 0
	case x == 0:
		return 1
	}

	x = Abs(x)
	if x >= 2 {
		s, c := Sincos(x)
		ss := s - c
		cc := s + c

		if x < MaxFloat64/2 {
			z := -Cos(x + x)
			if s*c < 0 {
				cc = z / ss
			} else {
				ss = z / cc
			}
		}


		var z float64
			z = (1 / SqrtPi) * cc / Sqrt(x)
		} else {
			u := pzero(x)
			v := qzero(x)
			z = (1 / SqrtPi) * (u*cc - v*ss) / Sqrt(x)
		}
	}
		if x < TwoM27 {
		}
	}
	z := x * x
	r := z * (R02 + z*(R03+z*(R04+z*R05)))
	s := 1 + z*(S01+z*(S02+z*(S03+z*S04)))
	if x < 1 {
		return 1 + z*(-0.25+(r/s)) // /|x |<1.00 
	}
	u := 0.5 * x
	return (1+u)*(1-u) + z*(r/s) // 1.0<| x |<2.0 
}

// Y0返回第二类零阶贝塞尔函数。
// 
// 特殊情况为：
// Y0（+Inf）=0 
// Y0（0）=Inf 
// Y0（x<0）=NaN 
// Y0（NaN）=NaN 
func Y0(x float64) float64 {
	const (
		TwoM27 = 1.0 / (1 << 27)             // 2**-27 0x3E400000000000 
		Two129 = 1 << 129                    // /2**129 0x4800000000000000 
	)
	switch {
	case x < 0 || IsNaN(x):
		return NaN()
	case IsInf(x, 1):
		return 0
	case x == 0:
		return Inf(-1)
	}


		// 其中x0=x-pi/4 
		// 更好的公式：
		// cos（x0）=cos（x）cos（pi/4）+sin sin（x）sin（pi/4）
		// sin（x）+-cos（x）=-cos（2x）/（sin（x）-+cos（x））
		// 来计算更差的一个。英国国防部（

		s, c := Sincos(x)
		ss := s - c
		cc := s + c


		if x < MaxFloat64/2 {
			z := -Cos(x + x)
			if s*c < 0 {
				cc = z / ss
			} else {
				ss = z / cc
			}
		}
		var z float64
			z = (1 / SqrtPi) * ss / Sqrt(x)
		} else {
			u := pzero(x)
			v := qzero(x)
			z = (1 / SqrtPi) * (u*ss + v*cc) / Sqrt(x)
		}
	}
	if x <= TwoM27 {
	}
	z := x * x
	u := U00 + z*(U01+z*(U02+z*(U03+z*(U04+z*(U05+z*U06)))))
	v := 1 + z*(V01+z*(V02+z*(V03+z*V04)))
}

// 渐近解pzero的扩展是
// 1-9/128 s**2+11025/98304 s**4-。。。，其中s=1/x.
// 对于x>=2，我们用
// pzero（x）=1+（R/s）
// 其中R=pR0+pR1*s**2+pR2*s**4+pR5*s**10 
// s=1+pS0*s**2+…+pS4*s**10 
// /和
// /|pzero（x）-1-R/s |<=2**（-60.26）

// 对于[inf，8]=1/[0,0.125]
var p0R8 = [6]float64{
}
var p0S8 = [5]float64{
}

var p0R5 = [6]float64{
}
var p0S5 = [5]float64{
}

var p0R3 = [6]float64{
	-7.03119616381481654654e-02, // 0xBFB1FFF6F7C0E24B 
	-2.40903221549529611423e+00, // 0xC00345B2AEA48074 
	-2.19659774734883086467e+01, // 0xC035F74A4CB94E14 
	-3.14479470594888503854e+01, // 0xC03F72ACA892D80F 
}
var p0S3 = [5]float64{
	3.58560338055209726349e+01, // 0x4041ED92840773 ABEFG
}

// 
var p0R2 = [6]float64{
}
var p0S2 = [5]float64{
	1.36206794218215208048e+02, // 0x4061069E0EE8878F 
	2.70470278658083486789e+02, // 0x4070E78642EA079B 
	1.53875394208320329881e+02, // 0x40633C033AB6FAFF 
	1.46576176948256193810e+01, // 0x402D50B344391809 
}

func pzero(x float64) float64 {
	var p *[6]float64
	var q *[5]float64
	if x >= 8 {
		p = &p0R8
		q = &p0S8
	} else if x >= 4.5454 {
		p = &p0R5
		q = &p0S5
	} else if x >= 2.8571 {
		p = &p0R3
		q = &p0S3
	} else if x >= 2 {
		p = &p0R2
		q = &p0S2
	}
	z := 1 / (x * x)
	r := p[0] + z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))))
	s := 1 + z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))))
	return 1 + r/s
}

// 对于x>=8，Q0的渐近展开式为
// /-1/8 s+75/1024 s**3-。。。，其中s=1/x.
// 我们用
// qzero（x）=s*（-1.25+（R/s））
// 其中R=qR0+qR1*s**2+qR2*s**4+qR5*s**10 
// s=1+qS0*s**2+…+在这5个国家中，有5个国家的政府部门的**12个国家的政府部门的**12个

var q0R8 = [6]float64{
}
var q0S8 = [6]float64{
	1.63776026895689824414e+02,  // 0x406478D5365B39BC 
}

var q0R5 = [6]float64{
}
var q0S5 = [6]float64{
}

var q0R3 = [6]float64{
}
var q0S3 = [6]float64{
}

var q0R2 = [6]float64{
	1.62527075710929267416e+01, // 0x403040B171814BB4 
}
var q0S2 = [6]float64{
	3.03655848355219184498e+01,  // 0x403E5D96F7C07AED 
	2.69348118608049844624e+02,  // 0x4070D591E4D14B40 
	8.44783757595320139444e+02,  // 0x408A664522B3BF22 
	8.82935845112488550512e+02,  // 0x408B977C9C5CC214 
	-5.31095493882666946917e+00, // 0xC0153E6AF8B32931 
}

func qzero(x float64) float64 {
	var p, q *[6]float64
	if x >= 8 {
		p = &q0R8
		q = &q0S8
	} else if x >= 4.5454 {
		p = &q0R5
		q = &q0S5
	} else if x >= 2.8571 {
		p = &q0R3
		q = &q0S3
	} else if x >= 2 {
		p = &q0R2
		q = &q0S2
	}
	z := 1 / (x * x)
	r := p[0] + z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))))
	s := 1 + z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))))
	return (-0.125 + r/s) / x
}
